Définir le sujet

Soit le sujet : Logements sociaux et qualification des chefs de ménages

Formuler des questions ou des hypothèses

Qu’elles soient justes ou fausses, les hypothèses permettent de cadrer l’analyse.

Organiser le travail

Sutout dans le cadre d’un groupe !

Charger les données statistiques

tab_ind<-readRDS("data/menag2018.RDS")

FALSE   COMMUNE   ARM      IRIS ACHL AEMM
FALSE 1   75056 75101 751010101 C114 2014
FALSE 2   75056 75101 751010101  A11 2008

Préparation de l’analyse

• Soit la relation entre logement en HLM (Y) et Diplôme le plus élevé du chef de ménage (X). Il s’agit de deux variables catégorielles (= qualitatives) que l’on va typiquement mettre en relation à l’aide d’un tableau de contingence et d’un test du chi-2. L’analyse statistique est simple sous R mais il faut tenir compte de trois difficultés

• Le choix de la population de référence est important. Ici on va sélectionner les ménages dont la personne de référence est âgée de 25-39 ans

• la sélection ou le regroupement des diplômes est également important car cela va influer sur les résultats du test.

• la pondération des individus doit également être prise en compte puisque le recensement est basé sur un sondage

Sélection des individus et des variables

#table(tab_ind$AGEMEN8)
tab_sel<- tab_ind %>% 
  filter(AGEMEN8 == "25") %>%
  select(DIPLM,HLML, IPONDL) 

Recodage des modalités

On cherche le code des modalités CS1 ezt HLML dans le fichier des métadonnées

meta<-readRDS("data/menag2018_meta.RDS")
metasel <- meta %>% filter(COD_VAR %in% c("DIPLM", "HLML"))
kable(metasel[,c(1,3,4)])

s

On recode les modalités des deux variables en regroupant certaines CSP

tab_sel$HLML<-as.factor(tab_sel$HLML)
levels(tab_sel$HLML)<-c("HLM-O","HLM-N",NA)
tab_sel$DIPLM<-as.factor(tab_sel$DIPLM)
levels(tab_sel$DIPLM) <- c("< BAC","< BAC","< BAC","< BAC","< BAC","< BAC",
                         "BAC","BAC",
                         "BAC+123","BAC+123","> BAC+3","> BAC+3",NA)
table(tab_sel$DIPLM)

FALSE 
FALSE   < BAC     BAC BAC+123 > BAC+3 
FALSE   65376   50352   88531  156494

Création du tableau de contingence non pondéré (FAUX)

La solution la plus simple semble être l’instruction table()

tab_cont<-table(tab_sel$HLML,tab_sel$DIPLM)

Création du tableau de contingence pondéré (JUSTE)

On pondère avec wtd.table() du package questionr.

library(questionr)
tab_cont_wtd<-wtd.table(tab_sel$HLML,tab_sel$DIPLM,
                        weights = tab_sel$IPONDL)

Comparaison des niveaux de dépendance automobile

• Tableau non pondéré … légèrement faux !

• Tableau pondéré … juste !

Visualisation du tableau de contingence

On choisit l’orientation du tableau et on l’affiche avec plot()

mytable<-wtd.table(tab_sel$DIPLM,tab_sel$HLML,weights = tab_sel$IPONDL)
plot(mytable)

Visualisation améliorée du tableau de contingence

Tant qu’à faire, on améliore la figure avec des paramètres supplémentaires :

Test du Chi-deux

Ce test se réalise facilement sur le tableau de contingence avec l’instruction chisq.test() :

mytest<-chisq.test(mytable)
mytest

    Pearson's Chi-squared test

data:  mytable
X-squared = 97192, df = 3, p-value < 2.2e-16

Visualisation des résidus

Lorsque la relation est significative, on visualise les cases les plus exceptionnelles avec mosaicplot( …, shade = T)

Le format sf (spatial features)

La cartographie et plus généralement les opérations géométriques sur des données spatiales dans R peuvent facilement être effectuées avec le package sf (spatial features) qui crée des objets uniques rassemblant à la fois

• un tableau de données (l’équivalent du fichier .dbf)
• une géométrie (l’équivalent du fichier .shp)
• une projection (l’équivalent du fichier .prj)

Lorsqu’on récupère des fonds de carte au format shapefile (.shp) ou dans d’autres formats standards comme GeoJson, la première tâche consiste donc à les convertir au formt sf afin de pouvoir les utiliser facilement dans R. L’importation se fait à l’aide de l’instruction st_read en indiquant juste le nom du fichier .shp à charger. Les autres fichiers (.dbf ou .proj) seront lus également et intégrés dans l’objet qui hérite de la double classe data.frame et sf.

Etapes de préparation des données

Dans notre exemple, nous allons suivre les étapes suivantes :

1. Préparer les données statistiques par IRIS dans un data.frame
2. Charger un fonds de carte par IRIS au format sf
3. Effectuer une jointure entre les deux fichiers par le code IRIS
4. Sauvegarder le résultat
5. Agréger les données statistiques et géométriques par commune
6. Sauvegarder le résultat.

Préparer les données statistiques

On importe le fichier des individus :

tab_ind<-readRDS("data/menag2018.RDS")

FALSE   COMMUNE   ARM      IRIS ACHL AEMM
FALSE 1   75056 75101 751010101 C114 2014
FALSE 2   75056 75101 751010101  A11 2008
FALSE 3   75056 75101 751010101  A11 2012

Agréger les données

On commence par créer un tableau long croisant les deux variables et leur effectif pondéré :

tab_long<- tab_ind %>%
           filter(HLML != "Y")%>%
           group_by(IRIS,HLML)%>%
           summarise(NB=sum(IPONDL))

Pivoter le tableau

Puis on fait “pivoter” le tableau pour l’obtenir en format large :

tab_large <- tab_long %>% pivot_wider(id_cols = IRIS, 
                                      names_from = HLML,
                                      names_prefix = "HLM_",
                                      values_from = NB,
                                      values_fill = 0)

Ajouter de nouvelles variables

On ajoute de nouvelles variables telles que le nombre total de ménage et le % de ménages en HLM :

tab<- tab_large %>% mutate(TOT = HLM_1+HLM_2,
                           HLM_pct = 100*HLM_1/TOT)

Examiner la distribution statistique

On examine l’histogramme donnant distribution statistique du % de ménages ordinaires résidant en HLM par IRIS.

p <- ggplot(tab) + aes (x = HLM_pct) +
                   geom_histogram(breaks = c(0,10,20,30,40,50,
                                             60,70,80,90, 100)) +
                   scale_x_continuous("% de ménages en HLM") +
                   scale_y_continuous("Nombre d'IRIS") +
                   ggtitle(label = "Distribution des logements sociaux dans Paris & PC",
                           subtitle = "Source : INSEE, RP 2018")

Charger les données géométriques

On importe le fichier des iris du Val-de-Marne qui est au format sf en ne gardant que les colonnes utiles

map_iris <- readRDS("data/map_iris.RDS")
map_iris<-map_iris[,c(4,5,1,2,7)]
names(map_iris)<-c("IRIS","NOM_IRIS","COM","NOM_COM","geometry")

FALSE [1] "sf"         "data.frame"

Visualisation du fonds iris avec sf

On peut facilement produire une carte vierge des iris du Grand Paris en faisant un plot de la colonne geometry du fichier sf

plot(map_iris$geometry,col="lightyellow")

Jointure des données IRIS et du fonds de carte

map_iris_tab<-merge(map_iris,tab,
                   by.x="IRIS",by.y="IRIS",
                   all.x=T,all.y=F)

Sauvegarde du fichier par IRIS

On sauvegarde notre fichier au format .RDS de R

saveRDS(map_iris_tab,"data/map_iris_hlm.RDS")

Agrégation statistique + géométriques

Grâce aux nouveaux packages de R (dplyr et sf) il est possible d’agréger simultanément les statistiques et les géométries après les avoir stockés dans un même objet de type “sf”

Du coup, on peut gagner beaucoup de temps dans les traitements et les analyses cartographiques, en particulier si l’on veut tester différents niveaux d’agrégation.

Agrégation des IRIS en communes

L’agrégation est très facile et elle concerne à la fois les variables (de stock) et les geometries

map_com_tab <- map_iris_tab %>% 
  group_by(COM, NOM_COM) %>% 
  summarise(HLM_1=sum(HLM_1,na.rm=T), 
            HLM_2=sum(HLM_2,na.rm=T)) %>%
  st_cast("MULTIPOLYGON")

map_com_tab <- map_com_tab %>%  mutate(TOT = HLM_1+HLM_2,
                                  HLM_pct = 100*HLM_1/TOT) 

Examiner la distribution statistique

On examine l’histogramme donnant distribution statistique du % de ménages ordinaires résidant en HLM par Commune.

p <- ggplot(map_com_tab) + aes (x = HLM_pct) +
                   geom_histogram(breaks = c(0,10,20,30,40,50,
                                             60,70,80,90,100)) +
                   scale_x_continuous("% de ménages en HLM") +
                   scale_y_continuous("Nombre de communes") +
                   ggtitle(label = "Distribution des logements sociaux dans Paris et PC",
                           subtitle = "Source : INSEE, RP 2018")

Sauvegarde du fichier par commune

On sauvegarde notre fichier au format .RDS de R

saveRDS(map_com_tab,"data/map_com_hlm.RDS")

Le package map_sf

Le package mapsf permet de réaliser des cartes statiques de très haute qualité. Il a en effet été mis au point par des cartographes et des géomaticiens professionnels de l’UMS RIATE. Il prend la suite du package cartography dont la maintenance demeurera assuré quelque temps encore mais ne fera plus l’objet de développements futurs. Le package mapsf présente l’avantage d’être totalement compatibvle avec le package sf ce qui n’était pas autant le cas pour le package cartography, plus ancien, et créé pour être compatible avec l’ancien package sp.

On trouvera la documentation du package mapsf à l’adresse suivante :

https://riatelab.github.io/mapsf/index.html

Création d’un template cartographique

Nous allons dans un premier temps apprendre à créer un fonds de carte vierge mais comportant tout l’habillage nécessaire (“template”). Pour cela nous allons charger différentes couches cartographiques correspondant respectivement au département, aux communes et aux iris :

map_iris<-readRDS("data/map_iris.RDS")
map_com <-readRDS("data/map_com.RDS")
map_dep <-readRDS("data/map_dep.RDS")

map_iris_hlm<-readRDS("data/map_iris_hlm.RDS")
map_com_hlm<-readRDS("data/map_com_hlm.RDS")

tracé d’un fonds de carte vierge

La fonction mf_map() avec le paramètre type = "base"permet de tracer une carte vide

 mf_map(map_iris, type = "base")

Superposition de couches

On peut toutefois ajouter toute une série de paramètres supplémentaire (col=, border=, lwd=, …) et superposer plusieurs fonds de carte avec le paramètre add = TRUE. L’ajout de la fonction layout permet de rajouter un cadre une légende.

# Trace les Iris avec des paramètres
mf_map(map_iris,  type = "base", 
       col = "lightyellow", border="gray50",lwd=0.3)
# Ajoute les contours des communes
mf_map(map_com,  type = "base", 
       col = NA,border="red",lwd=0.6,
       add = TRUE)
# Ajoute les contours des département
mf_map(map_dep,  type = "base", 
       col = NA,border="red",lwd=2,
       add = TRUE)
# Ajoute un cadre, un titre et des sources
mf_layout(title = "Paris & Petite Couronne", 
          credits = "Sources : IGN et INSEE")

Ajout d’un thème

On peut finalement modifier l’ensemble de la carte en lui ajoutant une instruction mf_theme() qui peut reprendre des styles existants ( “default”, “brutal”, “ink”, “dark”, “agolalight”, “candy”, “darkula”, “iceberg”, “green”, “nevermind”, “jsk”, “barcelona”) mais aussi créer ses propres thèmes

#Choix du thème
mf_theme("darkula")
# Trace les Iris avec des paramètres
mf_map(map_iris,  type = "base",
       border="white",
        lwd=0.3)
# Ajoute les contours des communes
mf_map(map_com,  type = "base", 
       col = NA, lwd=1,
       add = TRUE)
# Ajoute les contours des département
mf_map(map_dep,  type = "base", 
       col = NA, lwd=2,
       add = TRUE)
# Ajoute un cadre, un titre et des sources
mf_layout(title = "Paris & Petite Couronne", 
          credits = "Sources : IGN et INSEE")

Ajout de texte

On peut ajouter une couche de texte avec la fonction mf_label(). Par exemple, on va ajouter à la carte précédente le code insee des communes

mf_theme("agolalight")
mf_map(...)
mf_label(map_com, 
         var="NOM_COM",
         cex=0.4,
         col="blue",
         overlap = FALSE)
# Ajoute un cadre, un titre et des sources
mf_layout(title = "Communes et Iris du Val de Marne en 2017", 
          frame = TRUE,
          credits = "Sources : IGN et INSEE")

Carte de stock

Une carte de stock représente la localisation de quantités que l’on peut aditionner et dont le total a un sens. Par exemple un nombre d’habitants, un nombre de ménages, un nombre d’automobiles. Ce quantités doivent être représentées par des figures (cercles, carrés, …) dont la surface est proportionelle au stock afin que l’oeil du lecteur puisse les aditionner visuellement.

Dans le package mapsf, on réalise ce type de carte à l’aide de la fonction mf_map()en lui donnant le paramètre type="prop".

On va tenter à titre d’exemple de représenter la distribution du nombre de ménages ordinaires occupant un logement HLM par IRIS :

Carte de stock minimale

Les instructions minimales sont les suivantes :

# Trace les contours des communes
mf_map(x= map_iris, 
       type = "base")

# Ajoute le nombre de ménages par IRIS
mf_map(x =map_iris_hlm, 
      type ="prop",
      var = "HLM_1",
      add=TRUE)

Carte de stock habillée

mf_theme("agolalight")
mf_map(map_iris, type = "base",  
       col = "lightyellow",border="gray80", lwd=0.3)
mf_map(map_com, type = "base", 
       col = NA,border="black",lwd=1,add = TRUE)

mf_map(map_iris_hlm, var = "HLM_1",type = "prop",
  inches = 0.05, col = "red",leg_pos = "left",  
  leg_title = "Nombre de ménages", add=TRUE)

mf_layout(title = "Distribution des logements HLM en 2018", 
          frame = TRUE,
          credits = "Sources : IGN et INSEE")

Carte choroplèthe

Une carte choroplèthe ou d’intensité représente un phénomène relatif dont la somme n’a pas de sens. Par exemple, il serait absurde d’aditionner les % de logement HLM des IRIS du Val de Marne. Ces variables d’intensité caractèrisent donc l’état général d’une zone (choros) et elles vont être représentées par une couleur appliquée à toute la surface de la zone, d’où leur nom de cartes choroplèthes.

La fonction du package mapsf adaptée aux variables d’intensité est la fonction mf_map()munie du paramètre type = "choro".

On va prendre l’exemple du nombre de voitures par ménage.

Carte choroplèthe minimale

Si on ne précise rien, la carte est réalisée à l’aide de la palette par défaut avec un découpage des classes en quantiles (effectifs égaux).

# Carte choroplèthe
mf_map(
  x = map_iris_hlm, 
  var = "HLM_pct",
  type = "choro")

Carte choroplèthe habillée

On peut arriver à une carte beaucoup plus satisfaisante en contrôlant l’ensemble des paramètres de couleur et de découpage des classes. Puis en superposant les contours de communes au dessus de la carte des IRIS pour faciliter le repérage.

mybreaks = c(0, 10,20,30,40,50,
             60,70,80,90, 100)
mypal <- mf_get_pal(n = c(5, 5), 
                    pal = c("Greens", "Reds"))
# Carte choroplèthe des iris
mf_map( map_iris_hlm, var = "HLM_pct",
        type = "choro",
        breaks = mybreaks,pal = mypal, 
        border=NA,
       col_na = "gray80",leg_title = "% HLM",
       leg_val_rnd = 0)
# Contour des communes et cadre
mf_map(map_com, type = "base", col = NA,
       border="black",lwd=1,add = TRUE)
mf_layout(title = "% de ménages en HLM au RP  2018", 
          frame = TRUE,
          credits = "Sources : IGN et INSEE")

Carte stock + choroplèthe (1)

On peut combiner les deux modes cartographiques par superposition :

mf_theme("agolalight")

# Choisit les classes
mybreaks = c(0,5,10,20,40,80,100)

# Trace la carte choroplèthe
mf_map(
  x = map_iris_hlm, 
  var = "HLM_pct",
  breaks = mybreaks,
 # pal=mypal,
 type = "choro",
  border="white",
  col_na = "gray80",
 lwd=0.3,
 leg_title = "% ménages", 
 leg_val_rnd = 0,
  
)

# Ajoute les cercles proportionnels

mf_map(
  x =map_iris_hlm, 
  var = "HLM_1",
  type = "prop",
  inches = 0.06, 
  col = "red",
  leg_pos = "right",  
  leg_title = "Nb ménages",
  add=TRUE
)
# Ajoute les contours des communes
mf_map(map_com, 
       type = "base", 
       col = NA,
       border="black",
       lwd=1,
       add = TRUE)

# Ajoute un cadre, un titre et des sources
mf_layout(title = "Les ménages ordinaires en HLM  2018", 
          frame = TRUE,
          credits = "Sources : IGN et INSEE")

Carte stock + choroplèthe (2)

Mais les cercles dissimuent alors les plages de couleur, aussi on peut utiliser le type prop_choro qui place la variable choroplèthe à l’intérieur des cercles

mf_theme("agolalight")
mybreaks = c(0, 10,20,30,40,50,60,70,80,90, 100)
mypal <- mf_get_pal(n = c(5, 5), pal = c("Greens","Reds"))
mf_map(map_iris_hlm, type = "base",  
       col = "gray80",border="white", lwd=0.3)
mf_map(map_com, type = "base", 
       col = NA,border="white",lwd=1,add = TRUE)
mf_prop_choro( x = map_iris_hlm,  var = c("TOT", "HLM_pct"), 
  inches = 0.06, col_na = "grey", pal=mypal,
  breaks = mybreaks, nbreaks = 4, lwd = 0.1,
  leg_pos = c("right", "left"),leg_val_rnd = c(0,0),
  leg_title = c("nb. ménages", "% HLM"),
  add = TRUE)
mf_layout(title = "Les ménages ordinaires en HLM  2018", 
          frame = TRUE, credits = "Sources : IGN et INSEE")

Source

Le répertoire des logements locatifs des bailleurs sociaux (RPLS) a pour objectif de dresser l’état global du parc de logements locatifs de ces bailleurs sociaux au 1er janvier d’une année. Il est alimenté par les informations transmises par les bailleurs sociaux. La transmission des informations pour la mise à jour annuelle du répertoire des logements locatifs est obligatoire. Les données sont ensuite géolocalisées à l’adresse et mis à disposition des utilisateurs sur le site du ministère de la transition écologique

Les fichiers sont disponibles en général par régions mais livrés par départements dans le cas de l’Ile de France. Nous allons utilisé ici le fichier du 1er janvier 2020 accessible à l’adresse suivante

https://www.statistiques.developpement-durable.gouv.fr/le-parc-locatif-social-au-1er-janvier-2020-0

Métadonnées

Le fichier de données brutes au format .csv est accompagné d’un document excel précisant le code des variables et la façon dont elles ont été obtenues.

Spatialisation

Le fichier indique pour chaque logement sa localisation précise en terme d’adresse mais aussi d’étage dans un immeuble. A partir de ces données qualitatives, l’INSEE a procédé à un géocodage qui aboutit à la création de deux champs :

• coordonnées de latitude et longitude non projetées
• coordonnées de position en projection Lambert officielle

Selon les analyses on peut utiliser l’une ou l’autre de ces coordonnées. Mais la meilleur solution consiste à créer un fichier de type sf (spatial features) en coordonnées WGS94 qu’on pourra ensuite reprojeter dans le système de son choix.

Stratégie

Avant toute exploitation du fichier il est fortement recommandé d’analyser en détail les métadonnées et de définir une stratégie d’analyse.

1. choisir une première zone d’étude de petite taille et localisée de préférence dans un espace que l’on connaît bien.
2. choisir des variables intéressantes dont l’on connaît bien la signification et dont on a analysé en détail les métadonnées
3. vérifier la qualité des données en regardant notamment le nombre de valeurs manquantes, le dégré de précision, etc.
4. sélectionner des données auxiliaires issues d’autres sources que l’on souhaite croiser avec celles du RPLS en s’assurant de leur compatibilité (espace, temps, définition, …)
5. Ajouter les coordonnées spatiales et stocker le résultat dans un fichier de type sf comportant les indications de projection.

Importation du fichier

On importe le fichier enregistré au format RDS et on vérifie sa taille avec dim() et sont ype avec class()

don <- readRDS("data/RPLS2020.RDS")
dim(don)

[1] 844302     73

Le tableau comporte 844302 lignes (chacune correspondant à un logement) et 73 variables (décrites dans les métadonnées).

Choix de la zone d’étude

On décide de limiter notre analyse dans un premier temps à quatre communes voisines présentant des profils différents. On peut tester leur profil de respect de la loi SRU en 2019 sur l’application suivante : https://www.ecologie.gouv.fr/sru/

• Bonneuil-sur-Marne (94011): large excédent ( > 25% , pas de pénalité)
• Chennevières-sur-Marne (94019) : léger déficit (22.76%, 58 k€)
• Sucy-en-Brie: déficit (94071) : (19.93% , 150k€ )
• Ormesson-sur Marne (94055) : très fort déficit (2.29%, 665 k€)

On relève leur code INSEE afin de pouvoir faciliter l’extraction des données.

Choix des variables

On va se limiter ici à un très petit nombre de variables

Extraction du fichier

On applique la double sélection des individus et des variables en nous servant des fonctions filter()et select()du package dplyr.on aboutit ici à un fichier de 8139 lignes et 11 variables.

sel <- don %>% 
  filter(DEPCOM %in% c("94011","94019",
                       "94071","94055")) %>%
  select(result_id, result_label,
         DEPCOM, LIBCOM, 
        latitude,longitude,X,Y,
        CONSTRUCT,SURFHAB,NBPIECE) 

dim(sel)

[1] 8139   11

Recodage et typage

Certaines variables doivent être recodées ou changées de type afin de faciliter leur exploitation ultérieure par R.

sel$DEPCOM <- as.character(sel$DEPCOM)
sel$LIBCOM <- as.factor(sel$LIBCOM)
sel$PLG_IRIS <- paste(sel$DEPCOM,sel$PLG_IRIS, sep = "")
sel$SURFHAB <- as.numeric(sel$SURFHAB)

Résumé rapide

On analyse rapidement les variables thématiques choisies

Sauvegarde du fichier

On sauvegarde le fichier obtenu au format .RDS afin de garder le formatage des variables :

saveRDS(sel,"data/sel_logt.RDS")

Retour sur sf

Nous revenons sur le package sf (spatial features) que nous avons déjà rencontré au moment de la création de cartes thématiques par IRIS ou communes à l’aide du package mapsf.

Ici le package sf va être utilisé pour cartographier non pas des zones mais des localisations ponctuelles. Il pourra être à nouveau couplé avec le logiciel de cartogaphie statique comme mapsf , afin par exemple de placer les localisations des logements sociaux au dessus du fonds de carte des IRIS ou communes.

Mais il pourra aussi servir de base à des cartographies dynamiques permettant de placer les points sur des réseaux de rue et plus généralement sur des “tuiles” cartographiques permettant d’effectur des zoom. On utilisera à cet effet d’autres packages comme leaflet ou sa version simplifiée mapview.

Données ponctuelles

Nous reprenons le fichier de localisation établi au chapitre précédent et nous ne conservons que 6 variables:

logt <- readRDS("data/sel_logt.RDS") %>%
        select(adresse=result_id,
               X,Y,
               date = CONSTRUCT)

Données IRIS

Nous chargeons par ailleurs le fichier des IRIS en ne gardant que la zone d’étude :

map_iris <- readRDS("data/map_iris.RDS") %>%
         filter(INSEE_COM %in% c("94011","94019",
                       "94071","94055"))

Agrégation par commune

Rappel : on peut agréger les géométries d’un fonds sf. Ici on va créer le fonds de carte des communes.

map_com <- map_iris %>% group_by(INSEE_COM,NOM_COM) %>%
                summarise() %>%
                st_cast("MULTIPOLYGON")

Vérification de la projection

Nous savons que les coordonnées X,Y du fichier logement sont projetées en EPS 2154. Mais quelle est la projection de notre fonds IRIS ? S’agit-il de la même ?

st_crs(map_iris)$proj4string

[1] "+proj=lcc +lat_0=46.5 +lon_0=3 +lat_1=49 +lat_2=44 +x_0=700000 +y_0=6600000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"

st_crs(2154)$proj4string

[1] "+proj=lcc +lat_0=46.5 +lon_0=3 +lat_1=49 +lat_2=44 +x_0=700000 +y_0=6600000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"

A priori il s’agit bien de la même de sorte que les coordonnées X,Y devraient bien se superposer sur le fonds IRIS

Test de superposition

par(mar=c(0,0,0,0))
#trace les iris
plot(map_iris$geometry, 
     col="lightyellow", border="gray70",
     lwd=0.2)
# trace les communes     
plot(map_com$geometry, 
     col=NA, lwd=1, add=T)
# ajoute les points
points(x=logt$X,
       y=logt$Y, 
       cex=0.2,
       col="red",
       pch = 16)

fichier des adresses

Nous allons maintenant établir un fichier de localisation des adresses en nous servant de l’identifiant unique fourni par l’INSEE.

adr <- logt %>% select(adresse,X,Y) %>% 
               filter(duplicated(adresse) == F) %>%
               filter(is.na(X) ==F,is.na(Y)==F)

On constate qu’il n’y a que 652 adresses différentes alors que notre fichier fait état de 8139 logements. Une adresse regroupe donc en moyenne plus de 10 logements (habitat collectif).

Transformation en fichier sf

La transformation de notre fichier initial au format sf est facile à réaliser avec la fonction st_as_sf() du package sf. Mais il faut prendre garde de bien préciser le système de projection si l’on veut pouvoir ensuite l’utiliser.

map_adr <- st_as_sf(adr, coords = c("X","Y"))
st_crs(map_adr)<- 2154
str(map_adr)

Classes 'sf', 'data.table' and 'data.frame':    612 obs. of  2 variables:
 $ adresse : chr  "94019_0037_00005" "94019_0037_00001" "94019_0037_00003" "94019_0037_00007" ...
 $ geometry:sfc_POINT of length 612; first list element:  'XY' num  666780 6855840
 - attr(*, "sf_column")= chr "geometry"
 - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA
  ..- attr(*, "names")= chr "adresse"

Agrégation des logements

Notre nouveau fichier sf permet désormais d’effectuer des jointures avec le fichier des logements sociaux. A titre d’exemple on peut désormais compter le nombre de logements par adresse et leur ancienneté moyenne.

logt_by_adr <- logt %>% 
               group_by(adresse) %>%
               summarise(nblog = n(),
                         datemoy = mean(date))

Jointure

On peut désormais effectuer la jointure entre les données agrégées par adresse et le fichier sf de localisation des adresses :

map_logt <- inner_join(logt_by_adr,map_adr) %>% st_as_sf()

Cartographie avec mapsf

On peut désormais utiliser les méthodes de cartographie déjà vues avec mapsf :

mf_theme("agolalight")
mybreaks = c(1900, 1950, 1960, 1970, 1980, 1990, 2000, 2010, 2020)
mypal=brewer.pal(n = 8,name = "Spectral")
mf_map(map_iris, type = "base",  
       col = "gray80",border="white", lwd=0.3)
mf_map(map_com, type = "base", 
       col = NA,border="black",lwd=1,add = TRUE)
mf_prop_choro( x = map_logt,  var = c("nblog", "datemoy"), 
  inches = 0.08, col_na = "grey", pal=mypal,
  breaks = mybreaks, nbreaks = 4, lwd = 0.1,
  leg_pos = c("right", "left"),leg_val_rnd = c(0,0),
  leg_title = c("nb. logements", "ancienneté"),
  add = TRUE)
mf_layout(title = "Les logements sociaux en 2020",
        frame = TRUE, credits = "Sources : IGN et RPLS")

Sauvegarde des fichiers cartographiques

On sauvegarde nos différents fichiers cartographiques au format sf relatifs à la zone d’étude.

saveRDS(map_com,"data/sel_map_com.RDS")
saveRDS(map_iris,"data/sel_map_iris.RDS")
saveRDS(map_logt,"data/sel_map_logt.RDS")

Statique ou dynamique ?

• Cartographie statique
  • production d’images fixes de qualité
  • respect strict des règles de la sémiologie graphique
  • choix libre d’une projection adaptée (e.g. EPSG 2154)
  • production de documents imprimés à finalité normative ou scientifiques
• Cartographie dynamique
  • production d’interfaces consultables dans un navigateur.
  • modification possible de l’échelle et de l’arrière-plan
  • projection imposée par les “tuiles” (EPSG 4326)
  • production de documents interactifs à finalité citoyenne ou exploratoire

Packages R de cartographie dynamique

• leaflet : la référence
  • Une librairie javascript non liée à un langage (R, Python, html, …)
  • Disponible dans R sous forme de package
  • Développement constant
• ggmap : l’empire contre attaque
  • des outils cartogaphiques utilisant la syntaxe de tidyverse
  • impose désormais un lien avec Google
• tmap : une solution hybride
  • permet de passer facilement du mode statique au mode dynamique
• mapview : l’équivalent de mapsf
  • mis au point par des développeurs allemands
  • facilite l’usage de leaflet
  • en progrès constant (mais instable)

Préparation des données

On charge les fichiers au format sf et on les transforme en projection WGS94 (EPSG=4326), condition indispensable pour ajouter des “tuiles” dynamiques lors des zoom.

map_com <- readRDS("data/sel_map_com.RDS") %>%
              st_transform(4326)
map_iris <- readRDS("data/sel_map_iris.RDS") %>%
              st_transform(4326)
map_logt <- readRDS("data/sel_map_logt.RDS") %>%
              st_transform(4326)

Carte par défaut

Mapview produit par défaut une carte dynamique du fichier sf.

mapview(map_logt)

Superposition de couches

On peut créer des couches et les aditionner avec ‘+’ :

m1 = mapview(map_com, zcol = "NOM_COM") 
m2 = mapview(map_logt)
m1+m2

Exemple complet

On va essayer de reproduire la carte statique faite avec mapsf

# Carte des communes
map1 <- mapview(map_com, lwd=1, legend= FALSE,
                alpha.regions = 0.1)
# Carte des iris
map2 <- mapview(map_iris,lwd = 0.3, label= "NOM_IRIS",
                legend= FALSE, alpha.regions = 0)
# Carte des logements
map3 <- mapview(map_logt,
                zcol = "datemoy",
                at = c(1900,1960, 1970,1980,
                       1990,2000,2010, 2021),
                col.regions = brewer.pal(8, "Spectral"),
                cex= "nblog")
map1+map2+map3